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ABSTRACT 


We  develop  a computational  procedure,  based  on  Taylor’s  series 
and  continued  fractions,  for  evaluating  Tricomi's  incomplete  gamma 
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function  -y  (a,x)  = — — - Jet  dt  and  the  complementary  incomplete 


Ha 
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gamma  function  T(a,x)  = J e *ta  1 dt,  both  in  the  region  x > 0, 

x 

- OC  < Q < OD  . 
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AN  EVALUATION  PROCEDURE  FOR  INCOMPLETE  GAMMA  FUNCTIONS 

Walter  Gautschi 

1.  Introduction  . The  incomplete  gamma  function,  and  its  complementary  function,  are 
usually  defined  by 

X ® 

(1.1)  v(a.  x)  = / e lta  ' dt,  r(a,  x)  = / e ta  dt  . 

0 x 

By  Euler's  integral  for  the  gamma  function, 

(1.2)  v(a,x)  + r(a,x)  = T(a)  . 

We  are  interested  in  computing  both  functions  for  arbitrary  x,  a in  the  half-plane 

V = {(x,  a)  : x > 0,  - «>  < a < «> } . 

The  function  r(a,x)  is  meaningful  everywhere  in  V,  except  along  the  negative  a-axis, 
where  it  becomes  infinite.  The  definition  of  y(a,  x)  is  less  satisfactory,  inasmuch  as  it 
requires  a > 0.  The  difficulty,  however,  is  easily  resolved  by  adopting  Tricomi's  version  l 11 1 
of  the  incomplete  gamma  function, 

* . x a . . 

(1.3)  Y (a,  x)  = r(a)  y(a,  x)  , 

which  can  be  continued  analytically  into  the  entire  (x,  a)-plane,  resulting  in  an  entire  function 
both  in  a and  x, 

* . e~XM(  1 . a + l;x)  _ M(a.a  + !;-x) 

(1-4)  V (a.x>  = r(a  + 1)  r(a  + 1) 

where  M(a,b;z)  = 1 + g j |t  + " ' is  Summer's  function.  Moreover,  \ (a,x) 

is  real-valued  for  a and  x both  real,  in  contrast  to  F(a,x),  which  becomes  complex  for 
negative  x. 

Our  objective,  then,  is  to  compute  y (a,  x)  and  l(a,x)  to  any  piescribed  accuracy  for 

arbitrary  x,a  in  V.  We  do  not  attempt  here  to  compute  y(a,x)  for  negative  x,  which 
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may  well  be  a more  difficult  (but,  fortunately,  less  important)  task.  We  accomplish  our  task 
by  selecting  one  of  the  two  functions  as  primary  function,  to  be  computed  first,  and  computing 
the  other  in  terms  of  the  primary  function  by  means  of 
(1.5)  r(a,x)  = r(a){  1 - x%  (a,x); 

or 


(1.6)  y (a,  x)  - x a{l  - } ■ 

If  y (a,x)  is  the  primary  function,  we  evaluate  it  either  by  Taylor's  series,  or  by  a little- 

known  continued  fraction  due  to  Perron.  For  r(a,x)  we  use  exclusively  the  classical 

continued  fraction  of  Legendre.  Careful  attention  is  required  when  a is  very  close  (or  equa 

to  a nonpositive  integer,  in  which  case  (1.5)  is  not  suitable  for  computation. 

An  evaluation  procedure  of  the  generality  attempted  here  is  likely  to  be  of  interest  in 

* 


many  diverse  areas  of  application.  Widely-used  special  cases  of  \ 
include  K.  Pearson's  form  of  the  incomplete  gamma  function  [8], 

(1.7)  l(u,  p)  = (u\'p  + 1)P+1  y (p  + 1,  uVp  + 1),  u > 0, 
the  x^"Probability  distribution  functions 

(1.8)  P(x2lv)  = (jxY7  Y*(7,  7X2),  Q(x2Iv)  = -7- 

2 2 2 nj) 


(a,x)  or 


p > -1  , 


r(a,  x) 


1) 


the  exponential  Integrals 

(1.9)  E (x)  = xv_1r(-v  + 1,  x) 

V 

(which,  for  v - -n,  a negative  integer,  yield  the  molecular  integrals  A^fx)  ( 5]),  and  the 
error  functions 

(1.10)  erf  x = xY  (^,  x2),  erfc  x = ^ T(^,  x2)  . 

* 

When  a is  integer-valued,  y (a,  x)  becomes^n  elementary  function, 

(1.11)  YVn,x)  = xn,  Y*(n  + l,x)  = x (r‘+,^il  - e Xen(x)),  n=0,l,2,..., 

v k 

where  e (x)  = 2j  x A*. 
n k = 0 


-2 


2.  Asymptotic  formulas.  It  is  useful  to  have  some  general  idea  as  to  the  magnitude  of 


the  (unctions  r(a,  x)  and  v*(a,x),  when  one  or  both  variables  are  large.  If  a is 
bounded,  and  x — «,  it  is  well  known  [ 12,  p.  17  4]  that 

(2.1)  r(a,  x)  ~ e”Xxd"i,  x - oo,  |a  I bounded. 

By  (1. 6),  this  implies 

(2.2)  y*(a,  x)  ~ x’a,  x - <*>,  la  I bounded  . 

Equally  simple  is  the  case  |x|  bounded  and  a — » (over  positive  values  of  a),  in  which 
case  [12,  p.  175] 

+ -x 

(2.3)  Y (a>x)  ~ r^TT)  ’ a^°°’  ^ bounded  > 

hence,  by  (1.  5), 

(2.4)  r(a,x)~r(a),  a - <»,  lx  I bounded. 

An  indication  of  the  behavior,  when  both  variables  are  large,  can  be  gained  by  setting 
x - pa,  p > 0 fixed,  and  letting  a - <®.  Laplace's  method,  applied  to  the  integrals  in  ll.'.K 
then  gives 

(2.5)  r(a,pa)~ 


and 


(2.6)  Y (a,  pa)  ~ 
Similarly, 

(2.7)  r(-a,  pa)  ~ 


r-(a),  0 < p < 1 , 

\ Ha),  P = 1, 


a - oo  , 


taeiV!! 

(1  + p)a 


, P > 1 » 


r e~ap 


(I  - p)Ha  + 1) 


, 0 < P < 1 , 


1 _a  i 

2 a , p - l, 


a - oo 


(Pa)'a,  p > 1 , 


. .-a  -ap 

Lac) £ 

(1  + p)a 


, 0 < p 
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from  which,  by  (l.t>)  and  the  reflection  formula  for  the  gamma  function, 


(2.8) 


* 

Y (-a, pa)  - 


(TTTk  r(a)e  lf  ^ 


Uniform  asymptotic  expansions  are  derived  in  ( 10). 


a * 0 (mod  1), 


-4 


3.  Choice  of  primary  function  . We  recall  from  (1.  5),  (1.  6)  that  either  of  the  two  functions 
r(a,x),  y (a,x)  can  be  expressed  in  terms  of  the  other  by  means  of 

(3.1)  Ha,x)  = r(a){l  - x\  (a,  x) },  y (a,*)  r x afl  - } • 

In  our  choice  of  primary  function,  we  are  guided  by  two  considerations:  the  numerical  stability 
of  (3.1),  and  computational  convenience. 

With  regard  to  numerical  stability,  we  must  be  careful  not  to  loose  excessively  in  accuracy 
when  we  perform  the  subtractions  indicated  in  curled  brackets  in  (3.1).  No  such  loss  occurs 
if  the  absolute  value  of  the  respective  difference  is  larger  than,  or  equal  to,  ~ . This  criterion 
Is  easily  expressed  in  terms  of  the  ratio 


(3.2) 


r(a,  x) 


r(a,x) 

Ha)  ' 


Indeed,  the  first  relation  in  (3.1)  is  stable  exactly  if  I r( a , x)  I > ~,  while  the  second  is  stable 

in  either  of  the  two  cases  r(a,  x)  > ~ and  r(a,  x)  < ■“ . As  a consequence,  an  ideal  choice 

£tl  the  primary  function  is.  y (a,  x)  H < r(a,  x)  < ~ , and  I (a,  x)  if  lr(a,  x)  I < ~ ; in  all 

remaining  cases  either  choice  is  satis  factory. 

For  the  practical  implementation  of  this  criterion,  consider  first  a > 0,  x > 0.  In  this 

case,  0<r(a,x)<l,  and  r(a,  x)  increases  monotonically  in  the  variable  a ([11,  p-  27b]). 

Since  lim  r(a,  x)  = 0 and,  by  (2.  4),  lim  r(a,x)=l,  there  is  a unique  curve  a = p(x) 

a — o a — • 00 

In  the  first  quadrant  x > 0,  a > 0,  along  which  r(a,  x)  = ~,  and  r(a,x)  < ~ depending  on 

whether  a < p(x).  Since,  by  (2 . 5),  r(x,  x)  ~ j as  x - °°,  we  have  p(x)  = x for  x large. 

By  numerical  computation  it  is  found  that  in  fact  p(x)  = x for  all  positive  x,  the  value  of 

p(x)  consistently  being  slightly  larger  than  x.  Ideally,  a good  choice  for  the  primary  function 

* 

then  is  r(a,x)  if  0 < a < x,  and  y(a,x)  if  a>x>0. 

Direct  computation  of  I^ajX),  when  a and  x are  both  relatively  small,  however, 

is  difficult,  owing  to  the  sing'ular  behavior  of  D(a,x)  near  a = 0,  x = 0.  In  contrast, 
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♦ 

V (a,x)  is  easily  computable  from  its  Taylor  series,  when  x is  small  and  a is  arbitrary. 

These  practical  considerations  lead  us  to  choose  r(a,  x)  ^s  primary  function  only  if  x > 1.  5 

and  a < x,  while  y (a,  x)  is  taken  as  primary  function  in  ail  other  cases . The  choice  of 

the  breakpoint  x = 1.5  is  further  motivated  in  Section  5.2. 

An  analogous  discussion  in  the  case  a < 0,  x > 0 is  complicated  by  the  more  involved 

♦ 

behavior  of  r(a,x),  and  also  by  the  presence  of  lines  along  which  y (a,x)  vanishes. 

Computational  convenience,  nevertheless,  dictates  the  same  choices  of  primary  functions,  as 

made  above  in  the  case  a > 0,  x > 0.  Thus,  if  we  divide  the  half-plane  V into  three  regions 

* 

as  shown  in  Figure  3.1,  labelled  I,  II,  III,  we  use  y (a,x)  as  primary  function  in  regions 

1 and  II,  and  T(a,x)  in  region  III.  The  partition  of  Figure  3.1  also  reflects  different  choices 

* 

of  computing  methods:  Taylor's  series  for  y (a,x)  in  region  I,  Perron's  continued  fraction 
in  region  II,  and  Legendre's  continued  fraction  in  region  III. 

Since  practical  considerations  led  us  to  deviate  from  the  ideal  choice  of  primary  function, 
some  loss  of  accuracy  will  be  unavoidable.  The  problem  is  most  pronounced  in  region  1 

when  a is  close  to  a nonpositive  integer.  It  also  occurs  in  the  lower  portion  of  region  III 

* 

near  zero-curves  of  y (a,x).  The  simplest  remedy,  in  most  cases,  consists  in  computing 
the  primary  and  secondary  function  in  double  precision  (when  available).  An  exception  to 
this  arises  in  region  I,  when  a is  very  close  (or  equal!)  to  a nonpositive  integer,  in  which 
case  we  need  to  compute  the  appropriate  limit  function.  This  is  further  discussed  in  Section  5. 
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Taylor's  scries  for  y (a,x).  The  following  series  expansions  follow  immediately 

from  (1.4), 


(4.1) 


'’<-”*.1*7,  2 rf'r&r-*"  l 


T(a  + n + 1) 


n = 0 n = 0 

When  a > - ~~ , we  prefer  the  first  expansion,  since  it  terminates  after  the  first  term  when 
a = 0.  Although  the  series  is  alternating,  its  sum  is  always  positive  (in  fact,  larger  than 
e X)  when  a > -1,  as  is  seen  from  the  first  relation  in  (1.4). 

a < " 2>  we  let  (a)  denote  the  integer  closest  to  a [closest  to,  and  less  than, 
a,  in  case  of  a tie],  and  define  m = -{a)  - 1,  < - a - {a}.  Then 

a + 1 = -m  + < , 

where  m > 0 is  an  integer,  and  U I < •"  . From  the  second  series  in  (4.1)  [or,  alternatively, 

♦ * _x 

from  the  recurrence  relation  y (a  - 1,  x)  = xY  (a,x)  + e /r(a)  ] , we  then  find 

0 + —7  + 

« - m (c  - m)(t  - m + 1) 


*.  . m+1  * . 

V (a,x)  = x y (c,x)  + 


(4-  3) 


r(a  + 1) 


} . 


(*  - m)(<  - m + 1)  • • • (c  - 1) 

II  1 £ 

* I < ~ , we  can  use  the  first  expansion  in  (4.1)  to  compute  y («,x)  in  (4.3).  Also, 
when  m is  very  large,  we  need  to  compute  only  as  many  terms  in  the  sum  between  curled 
brackets  as  are  required  to  give  a meaningful  contribution  to  within  the  desired  accuracy. 

The  reciprocal  of  T(a  + 1),  appearing  in  (4.  3),  can  be  expressed  in  terms  of  m and 


• as 


(4.4) 


1 _ , ..m  < T( m t 1 - c ) 

r(a  + i)  = r(i  + « ) r(i  - «)  • 


By  virtue  of  m > 0,  I*  I < ~,  this  requires  an  evaluation  routine  for  r(x)  that  needs  to  be 
effective  only  in  the  range  x > ~ . It  should  be  noted,  however,  that  when  a is  very  close 
to  a negative  integer,  « will  be  small  and  inaccurate,  causing  the  evaluation  in  (4.4)  to  be 
Inaccurate . 
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5.  Evaluation  of  r(a,x)  when  a isnear  a nonpositive  integer.  Special  analysis  is 
required  for  Ha.x)  when  a is  very  close  to  a nonpositive  integer,  and  0<x<1.5.  The 
primary  function  then  is  y*<a,x>,  and  Equation  (1.  5)  expresses  r(a,x)  as  a product  of 
two  quantities,  one  very  large  in  absolute  value,  and  the  other  very  small  - the  difference  of 
two  almost  equal  numbers.  Computationally,  this  is  unsatisfactory  and  can  lead  to  gross 
inaccuracies.  If  m>0  is  an  integer,  however,  it  is  known  that 

*5,1)  r(-m  x)  x'mE  (x),  x > 0 

m-f  i ’ ’ 

Where  Em+l(x)  is  exponential  integral  of  order  m+1  (cf.  [ 1,  Eq.  5. 1. 4] ).  It  is  of 

Interest,  therefore,  to  discuss  under  what  conditions  the  right-hand  side  of  (5.1)  is  an  accept 

able  approximation  to  r(a,x),  when  a = -m,  and  how  to  compute  x'mE  (x). 

m+1 

5.1.  Limit  behavior  of  r(a,X)  as  a--m.  Consider 


-t  -m+  c 1 


r(-m  + i ,x)  J e t " *dt,  x > 0 , 


where  m>0  is  an  Integer,  and  |<  | is  small,  say,  |e  | < - . We  have 


H-m  + «,x)  = / e*t{t'm”1  + t’m"1(t*  -l)}dt  x‘mEm+1(x)+  J~e  ‘t  m I(t€  - l)dt  , 


so  that 


H-m  + € ,X)  - x'mE  (x)  = / e‘Vm'V  - l)dt 


Using  Taylor's  theorem,  we  have 


‘ 1 nt  , , 0«  / nt 

1 -e  -1+e  • «/nt,  0 0(c  ,t>,  0<e<lt 


giving 


(5.2) 


r(-m  -t  « ,x)  - x"mE  (x)  = « / e-Vm-Je0</ntintdt. 


-t  -m-1  (h/nt. 
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We  consider  first  the  case  0 < x < 1.  Since 


0«  / nt  - |«  | / nt  - |c  | 
e < e t if  0 < t < L , 


int  < t - 1,  if  t > 1 , 


we  obtain  for  the  integral  in  (5.  2), 


/ ---dt  </  e-tt-m'1t-  141  |/nt  Jdt  + / e-tt-m-1t,l|(t-1,dt 
x x ] 


<x  ‘ |/nx|x  mE  (x)  + f e"lt  m(t  - l)dt  . 


The  integral  on  the  far  right  equals  e'1  if  m = 0,  e’1  - E.(l)  if  m = 1 and  E (1)  - E (1) 

1 ’ m-1  m 


if  m > 2.  Using  the  inequalities  ([1,  Eq.  5.1.19]) 


-x  -x 

e r / e 

< E (X)  < 


xT7  < En(x)  - x +~n  - 1*  X>0-  n J>2.3 


we  find 


/ e 1 1 m(t  - l)dt  < a 


e-\  a Ji 
m ’ m \ 2 


1 if  m = 0 


if  m = 1 


— : if  m > 2 , 

m - 1 


and  therefore,  by  (5.2)  and  (5.3), 


r(-m  4 ( x)  - x*  E (x)  III  a e"1 

ail-  < I.  f 

i * W1  l *' 


Since  x k'm+j(x)  U(-ni,x)  clearly  decreases,  as  x Increases,  we  have,  again  with  the 
help  of  (5.  4), 

x'mEm+1,x)>Em+l(1)>^>  0 < x < 1 , 


and  thus,  finally, 


-10- 


(5.5) 


T(-m  + « ,x)  - x E^^lx) 


- m_ 

x E (x) 
m+l 


< |<  ] { x 1 |/nx|  + (m  + Z)a  } 0 < x < 1 

- m 


where 


(5.6) 


r 


(m  + Z)a  - \ 
m 1 


3 

2 

2(m  * Z) 

(in  4 l)(m  - 1) 


if  m = 0, 
if  m 1, 

if  m > 2 


In  the  case  x > 1,  we  estimate  the  integral  in  (5.  2)  by 


/ 


dt 


< f e *t  m *t  ( *1  nt  dt  < f e lt  m(t  - l)dt  , 


- 1 - m , 


from  which 


mr 


n-m  4 t ,X)  - x Em4l(x) 


x'mE  (x) 

m4l 


x E ,(x)  x 
m + 1 


e *t  m(t  - l)dt  . 


The  integral  on  the  right  evaluates  to  xe  if  m 0,  to  e - Ej(x)  if  m 1, 

x ^*E  ,(x)  - x ^*E  (x)  if  m > 2.  Using  once  more  the  inequalities  in  (5. 

m-1  m - 


C x(x  4 1)  if  m 0, 


(5.7) 


r(-m  4 t ,x)  - x‘mE  ,(x) 


m4l 


X’mEm4l(x) 


< l«  I ' < 


x (x  4 2)  if  m 1 , 
x 4 1 


x(x  4 m 4 1){ 


1 


x+m-2  x+m 


if 


If  we  restrict  attention  to  the  interval  1 < x < 1.  5,  the  factor  of  |i  | in  (5 
bounded,  respectively,  by  ^ = 3.75,  ^=3.15,  anc*  ~~  = 4.821...,  so 


(5.8) 


-mr 


U(-m  4 c ,x)  - x E^4  j(x) 


< 5 |c  | , 1 < X < 1.  5 . 


and  to 
1 , we  find 

m > 2 . 

7 I can  be 


. -11- 


Since  the  approximation  x E^+j(x)  will  only  be  used  if  |i  | is  quite  small,  we 

shall  assume,  somewhat  arbitrarily,  that  |<  | < .001.  Then  x"mE  ,(x)  approximates 

m4 1 

T(-m  4 < ,x)  correctly  to  d significant  decimal  places,  whenever 


lb 


I*  I v X 


.001 


1 nx  | 4 (m  4 2)o  } < ^ 10  , 0 < x < 1 

m — / 


(5.9) 


5 It  I < -^  10  d,  1 < x < 1.  5 . 


5-2.  Evaluation  of  x Em+j(x).  We  propose  to  use  the  expansion  ( 1,  Eq.  5.1.12] 

°°  i J,  , , h -m  m- 

*0,4  1)4  y -l-x>  ^ ) l 4 X—  y 

' m 4 1 -4  (r  4 l)(m4  14  r) ! / m - 

r=0  J sr0 


(5.10)  x'mE  .(x) 

m4l  m! 


m4l 

-l  nx  4 , 


-m  m-1  s 

-XI  m 


s ! (m  - s)  ’ 


where  the  last  term  is  understood  to  be  zero  if  m = 0,  and  where 


m 

4-(m  4 1)  = -y  4 l - 

k-1  K 

Y - • 5772.  . . being  Euler's  constant.  It  is  important  to  make  sure  that  little  cancellation 
take  place  among  the  various  terms  on  the  right  of  (5.10).  Empirical  computations  convinced 
us  that  no  more  than  approximately  one  decimal  digit  is  lost  due  to  cancellation,  if  we 
restrict  x to  the  interval  [0,1.5],  the  cancellation  being  noticeable  only  near  x 1.5 
when  m 0,1,2  or  3.  This,  in  fact,  is  one  of  the  motivations  for  placing  the  break  at 
x - 1.5  in  the  partition  of  Figure  3.1.  (Another  motivation  is  the  progressively  slowei 
convergence  of  Legendre's  continued  fraction  for  r(a,x)  as  x decreases.) 

Assuming,  then,  0 < x < 1.  5,  it  is  easily  verified  that  all  partial  sums  of  the  infinite 
series  in  (5.10),  as  well  as  the  sum  itself,  are  strictly  positive.  Since,  moreover, 

-fnx  4 ij,(m  4 1)>  -/n  1.5  + il<(2)  = - .4054...  4 .4227...  >0  for  m > 1,  we  see  that  the 
expression  between  curled  brackets  in  (5.10)  is  also  positive  for  m > 1.  No  such  statement 
can  be  made  for  the  finite  series  in  (5.10),  which  can  be  negative,  zero,  or  positive.  Its 
terms,  when  m is  large,  need  to  be  computed  only  so  long  as  they  contribute  significantly* 
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6.  Legendre's  continued  fraction  for  r<a,x).  The  following  continued  fraction,  due  to 
Legendre,  Is  well-known  ([  9 , p.  103],  [1,  Eg.  6.5.31]), 


(6.1) 


-a  Xr,,  , 1 1 - a 1 2 - a 2 

x e l (a,x  = — — — — — 

’ x+  1+  x+  1+ 


x+  1+  x+ 

It  converges  for  any  x > 0 and  for  arbitrary  real  a.  We  can  write  (6.1)  in  contracted  form  as 

3 3 3 

-a  xn,  , K0  K1  '2 

x 6 r(a>x)  rrrrr  tttt 
0 1 2 


0)c  = 2k  4 1 - a,  k = 0,1,2,...  , 

Pn  r 1,  Pw  = * k),  k = 1,2,3 


or,  alternatively,  in  the  form 


(6.2) 

where 

(6.3) 


-a  x 


(x  4 1 - a)x  e r(a,x)  = — — 


l a i !i  !i 

1+  14 


k(a  - k) 

ak  (x  4 2k  - 1 - a)(x  4 2k  4 1 - a) 


k - 1,2,3 


We  investigate  the  convergence  character  of  the  continued  fraction  in  (6. 2)  for  - <®  < a < x, 
x > 1.  5,  which  is  the  region  III  in  Figure  3.1,  where  (6 . 2)  is  going  to  be  used. 

It  is  well-known  (cf.  , e.g.,  [13,  p.  17 ff]  ) that  any  continued  fraction  of  the  form  (6 . 2) 
can  be  evaluated  as  an  infinite  series, 


(6.4) 

where 

(6.5) 


a.  a_  a. 


J_  _I  _2  _3  y 

14  14  14  14  ' Li 


k r 0 


P0  =1>  Pk  P1P2  pk>  k r1*2-3* 


(6.6) 


Pn  = °»  Pv  = 


~ak(M  pk-l’ 

1 + ak(1  + pk-l'  * 


k = 1,2,3,. 
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The  n-th  partial  sum  in  (6.4),  in  fact,  is  equal  to  the  n-th  convergent  of  the  continued  frac- 
tion, n 1,  2,  3, . . . . If  we  let 

% = 1 + pk  ’ 


then  the  recursion  for  p in  (6.6)  translates  into  the  following  recursion  for  o 


(6.7) 


° 0 *’  °k  1 4 a.o 

k k-1 


k 1,2,3,. 


Consider  now  the  case  of  a^  as  given  in  (6.  3).  If  k < a (thus  a > 1),  then  a^  > 0 
(since  a < x),  and  it  follows  inductively  from  (6. 7)  that  0 < < 1,  hence  -1  < p <•  0. 

In  view  of  (6.5),  this  means  that  (6.4)  initially  behaves  like  an  alternating  series  with  terms 

decreasing  monotonically  in  absolute  value. 

If  k > a,  then  a^  < 0,  and  may  become  larger  than  one.  However,  if 
0 < ofc  j < 2,  we  claim  that  1 < < 2 whenever  x > — . Indeed,  for  the  upper  bound  we 

must  show  that  1 4 a^a^  j > ~ , i.e.  , a^o^  j > ‘ 2 > or>  equivalently,  |a^  lo^  < . 

Since  j < 2,  it  suffices  to  show  la^l  < which  is  equivalent  to  1 < (x  - a)^  4 4kx. 
Since  k > 1 and  x > 0,  the  latter  is  certainly  true  if  x > ^ , which  proves  the  assertion 
<rk  < 2.  The  other  inequality,  1 < o^,  is  an  easy  consequence  of  1 4 t\ak  j > established 
in  the  course  of  the  argument  just  given,  and  the  negativity  of  a^  Since  for  the  first  k 
with  k > a we  have  0 < <r  ^ < 1 (by  virtue  of  the  discussion  in  the  preceding  paragraph, 
or  by  virtue  of  oQ  1),  it  follows  inductively  that  1 < < 2 for  all  k > a,  hence 

0 < p^  < 1.  In  the  case  k = a,  we  have  afc  = 0 and  = 1,  thus  pfc  0,  and  the 

argument  again  applies. 

We  have  shown  that  |pfc  I < 1 for  all  k > 1 , that  is,  the  terms  in  the  series  of  (6 . 4) 
are  nonincreasing  in  modulus,  whenever  -°°<a<x,  x>^,  in  particular,  therefore,  when 
(x,a)  is  In  the  region  III  of  Figure  3.1.  Moreover,  the  series  changes  from  an  alternating 
series,  initially,  to  a monotone  series,  ultimately. 
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In  the  region  a > x,  convergence  of  Legendre's  continued  fraction  may  deteriorate 
considerably  In  speed,  which,  together  with  the  appropriate  choice  of  primary  function,  is 
the  reason  why  we  prefer  Perron's  continued  fraction  for  a > x (cf.  Section  7). 
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7.  Perron's  continued  fraction  for  y (a,x).  From  the  differential  equation  satisfied  by 
Kummer's  function  M(a,b;z),  Perron  [9,  p.  278]  derives  the  following  continued  fraction, 


M(a  l.b  * 1 : z ) 
M(a,b;z) 


b (a  4 l)z  fa  ■»  2)z 
b - z4  b + 1 - z+  b 4 2 - z4 


We  have  shown  in  [ 2]  that,  when  z - -x  < 0 and  0 < a 4 1 < b,  the  continued  fraction  (7.1) 
converges  monotonically  in  the  sense  that  all  convergents  are  positive  and  increase  mono- 
tonlcally  to  their  limit.  Letting  b a + 1,  and  talcing  note  of  the  second  relation  in  (1.4), 
we  find  Immediately 


v (a  + 1.x) 

* 

Y (a,x) 


1 (a  4 lix  (a  < 2)x 

a + 1 4 x-  a+2-*x-  a+3  + x- 


❖ ♦ - x 

which  may  be  combined  with  xy  (a  + l,x)  y (a,x'  - e /Ha  + 1)  to  yield 

_.  r>,  . x , 1 x (a  + l)x  (a  + 2)x 

(7.2)  r(a  + De  y (a,x)  , v>  77TTT-  77^4  x-  a T3  V 7 

The  algorithm  in  (6.4)  - (6.6)  may  be  used  to  evaluate  <p , if  we  write 


a 4 1 4 x 


J_!i  !2  a_s 

1+  1+  1 4 1 4 


(a  4 k)x 

ak  (a  + k + x)(a  + k 4 1 4 x)  ’ 


k 1,2,3,...  . 


Perron's  continued  fraction  in  (7 . 2)  will  be  used  only  for  a > x > 1.  5,  i.  e.  , in  region  11  of 
Figure  3.1. 

The  speed  of  convergence  increases  with  a,  hence  is  slowest  when  a x.  Compared 

to  Legendre's  continued  fraction  (6.1),  Perron's  continued  fraction  unfortunately  converges 

only  about  half  as  fast  near  the  line  a x.  This  is  illustrated  in  Table  7.1,  which  shows 

the  number  of  convergents  required  for  relative  accuracies  ^10  d,  d = 8(4)20,  Legendre's 

continued  fraction  (L)  being  used  for  a - .999xr,  Perron's  continued  fraction  (P)  for 

a =1.001x  , where  x = . 3872  4 5(r-l),  r = 5(5)40. 
r’  r ’ 
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D 

5 

10 

15 

20 

25 

30 

35 

40 

1 

B 

20. 3872 

46. 3872 

70. 3872 

9 5.  387  2 

120. 3672 

146. 3872 

170. 3872 

195.  3672 

8 

L 

16 

20 

24 

27 

29 

31 

33 

34 

P 

29 

40 

49 

55 

61 

66 

71 

76 

12 

L 

18 

25 

30 

34 

37 

40 

42 

44 

P 

38 

53 

64 

73 

81 

88 

94 

100 

16 

L 

20 

30 

36 

40 

44 

47 

50 

52 

P 

46 

64 

77 

88 

97 

106 

113 

121 

20 

L 

22 

34 

41 

46 

50 

54 

57 

60 

P 1 

54 

74 

89 

101 

112 

122 

1 30 

1 39 

Table  7.1.  Speed  of  convergence  of  Legendre's  (L)  and  Perron's  (P)  continued 
fraction  near  the  line  a x. 


Equation  (7.2)  is  subject  to  a potential  loss  of  accuracy  when  <p  L 1 . From  the 

asymptotic  results  in  Section  2,  however,  we  can  see  that  <p  — 0 when  x is  fixed  and 

a -»  <*>.  When  x = pa  and  p is  fixed  in  0 < p < 1,  then  </’-*■  p as  a — 00 . Any 

significant  cancellation,  therefore,  appears  to  be  possible  only  near  the  line  x a,  when 

* 1 -a 

a is  large.  Even  then,  however,  the  problem  is  very  minor,  since  y (a, a)  ~ — a , hence 

1 - e - 1 -./—  (1  + o(l)).  as  a -»  «>.  We  thus  loose  between  one  and  two  decimal  digits 
V "a  ’ 

. * . -1800 

only  if  a ~ 640,  at  which  point  y (a, a)  = 10  , much  too  small  to  be  representable 

on  most  computers. 
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8.  Flowchart  . We  briefly  review  the  salient  features  of  the  complete  evaluation 


procedure,  which  forms  the  basis  of  the  algorithm  in  [ 3],  and  encapsulate  them  in  the  flowchart 
of  Figure  8.  1 below. 

* 

Our  goal  is  to  evaluate  y (a,  x)  and  I'(a,x)  to  d correct  significant  decimal  places, 
where  d < d^,  with  dQ  denoting  the  number  of  decimal  digits  available  on  the  computing  device 
at  hand.  (On  a binary  computer  with  (3  binary  digits  in  the  mantissa  of  the  single-precision 
floating-point  word,  we  may  set  dQ  = (3/n2//nlO.)  Due  to  unavoidable  rour.ding  errors, 
and  the  possibility  of  minor  cancellation  errors  at  a few  places  (e.g. , in  Taylor's  series 
for  the  exponential  integral  and  in  Perron's  continued  fraction;  cf.  Section  5.2  and  the  final 
paragraph  of  Section  7),  it  would  be  unrealistic  to  require  d = dQ . We  shall  in  fact  tolerate 
losses  of  1 to  3 decimal  digits,  before  taking  any  precautions,  so  that  a reasonable 
choice  of  d should  be  less  than  dQ  - 1 or  d - 2. 

Due  to  our  choice  of  primary  function  (cf.  Section  3),  more  serious  cancellation  errors 
may  be  incurred  in  the  use  of 


(8-1) 

when  0 < x < 1.  5 and 
(8.2) 


T(a,x)  = F(a){l  - xaY  (a,  x) } 
a < x,  and  in  the  use  of 


Y% 


when  x > 1.5  and  a < 0.  With  regard  to  (8.2),  the  problem  is  only  significant  near  a zero 
♦ 

of  y (a,x).  In  our  algorithm  in  ( 3]  we  make  no  attempt  at  achieving  full  significant  accuracy 
♦ 

near  zero-curves  of  y (a,x).  Instead,  we  issue  an  appropriate  message  to  the  user  of  the 
algorithm,  if  the  accuracy  requirement  cannot  be  met. 

In  (8.1),  the  loss  of  accuracy  is  most  prominent  when  a is  close  to  a nonpositive  integer, 

a = -m,  m = 0,  l,  2,  . . . . 

Different  actions  are  taken,  in  this  instance,  depending  on  whether  double-precision  Is  avail- 
able or  not.  If  it  is,  we  recompute  the  expression  on  the  right  of  (8.  1)  in  double  precision, 
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whenever  the  single-precision  calculation  (to  d digits!)  reveals  that 


(«•  3)  ll  - xaY  (a,  x)  I < .05  , 

1 . e. , when  the  d-digit  accuracy  in  I(a,x)  is  in  jeopardy.  If  single  precision  is  the  only 

d-d 

precision  mode  available,  and  if  enough  extra  precision  is  available  (say,  10  > 20),  the 

occurrence  of  (a.  3)  will  prompt  us  to  recompute  I (a,  x)  to  full  single-precision  accuracy. 

If  these  countermeasures  prove  to  be  of  no  avail,  then  again  our  algorithm  issues  an  appropriate 
message.  This  will  notably  be  the  case  if 

= * -(d-d) 


(*•4) 


11  - x y (a, x)  I < 10 


after  y has  been  recomputed  to  an  accuracy  of  dQ  decimal  digits. 

It  is  possible  to  estimate  how  close  to  a nonpositive  integer  -m  the  value  of  a has 

to  be  for  the  predicament  (a.  4)  to  occur.  From  (8.  2),  and  r(-m,x)  = x"mE  ,(x),  it  follows 

m+1 

Indeed  that 

V (a,  x)  = x - [x  tnx  + (-1)  m!Em+1(*)]<  + 0(c  ),  e = a + m - 0 , 


from  which 


1 - x“y  (a,  x)  = (-l)mm!x  mE  (x)c  + 0(tc),  c - 0 . 

m+i 


Neglecting  the  0(<  )-term,  we  thus  have  (a . 4)  if 


„ -(d-d) 

m!x  Em+J(x)  U I < 10 


which,  by  virtue  of  x~  Em+,(x)  > 1.  5_mEmi](  1 . 5),  implies 
(«.5)  |<  | < k 10  (°° 


m + l 
-(dn-d) 

1( 

m m 


m!1-5  Em  + l<‘-5) 


, 226. 


(The  first  few  values  of  are:  k = 9 .9980,  k = 20 . 520,  = 19 . 827 , k =12. 

m U 1 2 3 

k4  = 5.  4746,  = 1.9137,...  .)  Thus,  for  our  procedure  to  concede  defeat,  the  value  of  a 

must  be  sufficiently  close  to  -m  so  as  to  have  (a.  5),  but  not  so  close  as  to  satisfy  (5.  9), 
since  then  the  limit  value  (for  a * -m)  will  be  accurate  enough. 
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figure#.!.  Flow  chart  for  computing  y*(aix)  and  Ha,  x) 
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An  important  feature  of  our  algorithm  is  the  automatic  monitoring  of  overflow  and  under- 
flow conditions.  This  is  accomplished  by  first  computing  the  logarithm  of  the  desired 
quantities  and  by  making  the  tests  for  overflow  and  underflow  on  the  logarithms.  As  a result, 
minor  inaccuracies  are  introduced  in  the  final  exponentiation,  particularly  if  the  result  is 
near  the  overflow  or  underflow  limit. 
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9.  Testing • The  algorithm  in  [ 3 J , and  a double-precision  version  of  it,  were  tested 
extensively  on  the  CDC  6S00  computer  at  Purdue  University  and  on  the  UNIVAC  1110  computer 
at  the  University  of  Wisconsin.  The  double-precision  algorithm  was  used  to  provide  reference 
values  for  checking  the  single-precision  algorithm,  and  on  a few  occasions,  to  check  against 
high-precision  tables  (notably  the  MS  tables  in  f 6]).  Other  reference  values  were  taken 
from  various  mathematical  tables  in  the  literature. 

The  tests  include: 

(i)  the  error  functions  (1  . 1 0),  checked  against  Tables  7 . 1 and  7 . 3 in  [ 1 J ; 

(ii)  the  case  ( 1 . 1 1 ) of  integer  values  a = n,  -20  < n < 20; 

(iii)  the  exponential  integral  ( x)  in  (1.9)  for  integer  values  v = n,  0 < n < 20, 

and  fractional  values  of  v in  0 < v < 1,  checked  against  Tables  1,  II,  111  in  [ 7 ]; 

(iv)  Pearson's  incomplete  gamma  function  (1.7),  checked  against  Tables  I and  11  in  [ej; 

(v)  the  incomplete  gamma  function  P(a,x)  = (x/2)a  y*(a,  x/2),  checked  against  the 

tables  in  (4); 

(vi)  the  x2-distribution  (1. 8),  checked  against  Table  26. 7 in  [ 1 J ; 

(vii)  the  molecular  integral  An(x),  checked  against  Table  1 in  ( S)  and  the  more  accurate 
tables  in  ( 6) . 

In  addition,  the  performance  of  the  routines  was  examined  for  parameter  values  a in 
the  neighborhood  of  zero  and  in  the  neighborhood  of  negative  integers.  When  a is  very  close 
to  a negative  integer,  the  accuracy  of  y (a,x)  was  observed  to  deteriorate  substantially. 

This  is  due  to  inaccuracies  in  the  evaluation  of  T(a  + 1),  as  remarked  at  the  end  of  Section  4. 
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An  important  feature  of  our  algorithm  is  the  automatic  monitoring  of  overflow  and  under- 
flow conditions.  This  is  accomplished  by  first  computing  the  logarithm  of  the  desired 
quantities  and  by  making  the  tests  for  overflow  and  underflow  on  the  logarithms.  As  a result, 
minor  inaccuracies  are  introduced  in  the  final  exponentiation,  particularly  if  the  result  is 
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